Introduction

The first official case of COVID-19 in the USA has been confirmed on the 21st of January. About three months later, almost 1 million cases have been discovered. In this context of this pandemic, it is of utmost importance to understand how the pandemic evolves by reporting data in a clear an insgihtful way.

This project has two main goals:

visualize different COVID-related metrics: infection rate / 100k individuals, total cases, and total deaths at the county level (by representing maps and time series)
identify counties with potential errors in official counts: it has been shown that some counties negative counts in cumulative cases, which is not possible.

Notes about the report:

R code used to generate this report is provided. You just need to click on the “Code” button on the right to display the code used in a given section.
Most graphs and tables are interactive: you can zoom in and out, click on elements to display more content, or search for specific data points.

Data preparation

Dependencies and input files

We first need to load R packages (mostly used for interactive visualizations).

Then we set the location of input files used in this project:

a web scraped data file (april15.csv) containing COVID data for a single time-point
official population counts in US counties (census), directly taken from the web
the New-York Times COVID-19 data, taken from the web
official counties borders to be displayed on the map (downloaded from census.gov)

## R packages to have installed
library(rgdal)
library(dplyr)
library(leaflet)
library(scales)
library(dygraphs)
library(crosstalk)
library(plotly)

## Local paths or URLs to files

# web-scraped data
CSV.input <- "./data/april15.csv"

# Census data
CENSUS.input <- "https://www2.census.gov/programs-surveys/popest/datasets/2010-2019/counties/totals/co-est2019-alldata.csv"

# NYT data
NYT.url <- "https://raw.githubusercontent.com/nytimes/covid-19-data/master/us-counties.csv"

# Counties borders
# downloaded from https://www.census.gov/geographies/mapping-files/time-series/geo/carto-boundary-file.html
# and put in a 'shp' subdirectory
MAP.dir <- "./shp"
MAP.prefix <- "cb_2018_us_county_20m"

Loading and cleaning the web scraped file

We first load the april15.csv file. We keep the unique county ID (FIPS or GEOID here), the number of cases and deaths.

### Begin data prep
dat <- read.csv(CSV.input, header = TRUE, sep=",", stringsAsFactors = FALSE)

# adding leading zeros to any FIPS code 
# that's less than 5 digits long to get a good match
dat <- data.frame(
  GEOID = formatC(dat$fips, width = 5, format = "d", flag = "0"),
  cases = dat$cases,
  deaths = dat$deaths
)

head(dat) # check how it looks

Loading and cleaning census data

Census data taken from census.gov will be useful to scale the number of cases with population estimates for each county. After loading the file, we generate a table with the properly formatted unique county ID (GEOID) and the 2019 population estimate.

### Get population data
census <- read.csv(CENSUS.input, sep = ",", header = TRUE)

# Prepare county IDs (concatenate state and county IDs)
census$GEOID <- paste0(
  formatC(census$STATE, width = 2, format = "d", flag = "0"),
  formatC(census$COUNTY, width = 3, format = "d", flag = "0")
)

# we keep only the county level from this dataset (encoded as SUMLEV = 50)
census <- census[census$SUMLEV == 50, ]

# we just keep the GEOID (= FIPS) and latest census estimate (2019)
census.small <- data.frame(
  GEOID = census$GEOID,
  pop = census$POPESTIMATE2019
)
rownames(census.small) <- census.small$GEOID

Getting and preparing NYT data

Another dataset we’re gonna look at is provided by the New-York Times. We load it as a data frame containing different variables: the date, county, state, FIPS (county ID), cases and deaths. We format the county ID just like before.

## Load data
nyt.data <- read.csv(
  NYT.url,
  stringsAsFactors = FALSE
)

## Format county ID
nyt.data$GEOID <- formatC(nyt.data$fips, width = 5, format = "d", flag = "0")

Compute rate per 100k using census data for both datasets

For the two COVID datasets, we have the number of cases and deaths. We can compute for each dataset two other metrics: the rate of cases per 100k inhabitants, and the rate of deaths per 100k.

dat$pop <- census.small[dat$GEOID, ]$pop
dat$cases100k <- 100000 * dat$cases / dat$pop
dat$deaths100k <- 100000 * dat$deaths / dat$pop

nyt.data$pop <- census.small[nyt.data$GEOID, ]$pop
nyt.data$cases100k <- 100000 * nyt.data$cases / nyt.data$pop
nyt.data$deaths100k <- 100000 * nyt.data$deaths / nyt.data$pop

Preparing the base map

We first load county borders and remove Alaska and islands.

## Load the shape file
us.map <- readOGR(dsn = MAP.dir, layer = MAP.prefix, stringsAsFactors = FALSE)

## IMPORTANT preprocessing step

# Remove Alaska(2), Hawaii(15), Puerto Rico (72), Guam (66), 
# Virgin Islands (78), American Samoa (60), Mariana Islands (69), 
# Micronesia (64), Marshall Islands (68), Palau (70), Minor Islands (74)
us.map <- us.map[!us.map$STATEFP %in% c("02", "15", "72", "66", "78", "60", "69",
                                        "64", "68", "70", "74"),]

# Make sure other islands are removed
us.map <- us.map[!us.map$STATEFP %in% c("81", "84", "86", "87", "89", "71", "76",
                                        "95", "79"),]

Maps based on web-scraped and NYT data

We want to represent our COVID data on a map of the USA. To do so, we will generate maps using the Leaflet framework. The idea is to add polygons representing counties to the basemap. Polygon coloring depends on the metric of interest (here, total cases or rate /100k). Clicking on a county gives more information about this area.

# Merge spatial df with downloade ddata.
leafmap <- merge(us.map, dat, by=c("GEOID"))

# Format popup data for leaflet map.
popup_dat <- paste0("<strong>County: </strong>", 
                    leafmap$NAME, 
                    "<br><strong>Cases: </strong>", 
                    round(leafmap$cases),
                    "<br><strong>Cases /100k: </strong>", 
                    round(leafmap$cases100k),
                    "<br><strong>Deaths: </strong>", 
                    round(leafmap$deaths),
                    "<br><strong>Deaths /100k: </strong>", 
                    round(leafmap$deaths100k))

save(leafmap, file = "./data/leafmap.rda")

Total number of cases (April 15 data)

First, we represent the total number of cases in each county.

pal <- colorQuantile("YlOrRd", leafmap$cases100k, n = 10, na.color = "#fff8e8")

### Render map in leaflet
leaflet(data = leafmap) %>% addTiles() %>%
  addPolygons(fillColor = ~pal(cases), 
              fillOpacity = 0.8, 
              color = "#BDBDC3", 
              weight = 1,
              popup = popup_dat)

Rate per 100,000 (April 15 data)

Second, we represent the rate of cases per 100,000 individuals.

pal <- colorQuantile("YlOrRd", leafmap$cases100k, n = 10, na.color = "#fff8e8")

### Render map in leaflet
leaflet(data = leafmap) %>% addTiles() %>%
  addPolygons(fillColor = ~pal(cases100k), 
              fillOpacity = 0.8, 
              color = "#BDBDC3", 
              weight = 1,
              popup = popup_dat)

Prevalence with (NYT data, last day available)

# We first create a data frame containing data from the last day
last.day <- (nyt.data[nyt.data$date == max(as.Date(nyt.data$date)), ])
last.day <- last.day[order(last.day$cases, decreasing = TRUE), ]

# Merge spatial df with downloade ddata.
leaf.nyt <- merge(us.map, last.day, by=c("GEOID"))

# Format popup data for leaflet map.
popup_dat <- paste0("<strong>County: </strong>", 
                    leaf.nyt$NAME, 
                    "<br><strong>Cases: </strong>", 
                    round(leaf.nyt$cases),
                    "<br><strong>Cases /100k: </strong>", 
                    round(leaf.nyt$cases100k),
                    "<br><strong>Deaths: </strong>", 
                    round(leaf.nyt$deaths),
                    "<br><strong>Deaths /100k: </strong>", 
                    round(leaf.nyt$deaths100k))

save(leaf.nyt, file = "./data/leafnyt.rda")
pal <- colorQuantile("YlOrRd", leaf.nyt$cases100k, n = 10, na.color = "#fff8e8")

### Render map in leaflet
leaflet(data = leaf.nyt) %>% addTiles() %>%
  addPolygons(fillColor = ~pal(cases100k), 
              fillOpacity = 0.8, 
              color = "#BDBDC3", 
              weight = 1,
              popup = popup_dat)

## Format county ID
fips.20 <- last.day$GEOID[2:21]
fips.100 <- last.day$GEOID[2:101]

Time series representation based on the New York Times data

We can now have a look at COVID data over time. To do so, we will represent COVID cases evolution for the 100 most affected counties (most affected = highest number of cases so far).

Total cases over time

sd <- SharedData$new(nyt.data[nyt.data$GEOID %in% fips.100,], key = ~county)
fig <- plot_ly(sd, x = ~(date)) 
fig <- fig %>% add_trace(
  y = ~cases,
  color = ~GEOID,
  name= ~county,
  type = 'scatter',
  mode = 'lines+markers',
  marker = list(size = 2),
  line = list(width = 1))  

bscols(
  list(
    # filter_slider("mag", "Filter cases", sd, column = ~cases, step=2500, width=250),
    filter_select("stt", "Select state name:", sd, ~state),
    filter_select("ctw", "Select county name:", sd, ~county),
    filter_select("fps", "Select county FIPS:", sd, ~GEOID)
  ),
 fig,
  widths = c(3, NA)
)

Select state name:

Select county name:

Select county FIPS:

Rate of infection over time

sd <- SharedData$new(nyt.data[nyt.data$GEOID %in% fips.100,], key = ~county)
fig <- plot_ly(sd, x = ~(date)) 
fig <- fig %>% add_trace(
  y = ~cases100k,
  color = ~GEOID,
  name= ~county,
  type = 'scatter',
  mode = 'lines+markers',
  marker = list(size = 2),
  line = list(width = 1))  

bscols(
  list(
    # filter_slider("mag", "Filter cases", sd, column = ~cases, step=2500, width=250),
    filter_select("stt", "Select state name:", sd, ~state),
    filter_select("ctw", "Select county name:", sd, ~county),
    filter_select("fps", "Select county FIPS:", sd, ~GEOID)
  ),
 fig,
  widths = c(3, NA)
)

Select state name:

Select county name:

Select county FIPS:

Identifying counties with negative numbers

Data table

nyt.data.sorted <- nyt.data[order(nyt.data$date), ]
diff.per.county <- tapply(nyt.data.sorted$cases, nyt.data.sorted$GEOID, diff)
diff.per.county.max <- lapply(diff.per.county, min)
differences <- data.frame(GEOID = names(diff.per.county.max), Max.diff = unlist(diff.per.county.max), stringsAsFactors = FALSE)
differences <- differences[order(differences$Max.diff), ]

## Add county and state

DT::datatable(differences[differences$Max.diff < 0 & differences$GEOID != "   NA", ])

Time series for counties with largest discrepancies

which.diff <- lapply(diff.per.county, function(x) any(x < (-10)))
which.diff <- do.call(rbind, which.diff)
names.diff <- names(which.diff[which.diff[, 1] == TRUE, ])[-1]

fips.diff <- leafmap$GEOID[leafmap$GEOID %in% names.diff]

print(paste(leafmap$NAME[leafmap$GEOID %in% names.diff], fips.diff))

##  [1] "Cullman 01043"     "Onondaga 36067"    "Tazewell 17179"   
##  [4] "Dougherty 13095"   "Carson City 32510" "Ripley 18137"     
##  [7] "Oakland 26125"     "Madison 01089"     "Lafayette 22055"  
## [10] "Rensselaer 36083"  "St. Charles 22089" "Tuscaloosa 01125" 
## [13] "Lexington 45063"   "St. Landry 22097"

sd <- SharedData$new(nyt.data[nyt.data$GEOID %in% fips.diff,], key = ~county)
fig <- plot_ly(sd, x = ~(date)) 
fig <- fig %>% add_trace(
  y = ~cases,
  color = ~GEOID,
  name= ~county,
  type = 'scatter',
  mode = 'lines+markers',
  marker = list(size = 2),
  line = list(width = 1))  

bscols(
  list(
    # filter_slider("mag", "Filter cases", sd, column = ~cases, step=2500, width=250),
    filter_select("stt", "Select state name:", sd, ~state),
    filter_select("ctw", "Select county name:", sd, ~county),
    filter_select("fps", "Select county FIPS:", sd, ~GEOID)
  ),
 fig,
  widths = c(3, NA)
)

Select state name:

Select county name:

Select county FIPS:

Represent counties with discrepancies on a map

Work in progress.

COVID-19 in US counties

COVID-19 in US counties

Introduction

Data preparation

Dependencies and input files

Loading and cleaning the web scraped file

Loading and cleaning census data

Getting and preparing NYT data

Compute rate per 100k using census data for both datasets

Preparing the base map

Maps based on web-scraped and NYT data

Total number of cases (April 15 data)

Rate per 100,000 (April 15 data)

Prevalence with (NYT data, last day available)

Time series representation based on the New York Times data

Total cases over time

Rate of infection over time

Identifying counties with negative numbers

Data table

Time series for counties with largest discrepancies

Represent counties with discrepancies on a map